Background:
Author: Katherine Wynne (Wynnekat@msu.edu)
Co-authors: E. Eyerly, L. Sullivan
Created: 14 June 2023
Files:
Cleaned_Spring_Cover_Data_July_2019.xlsx - contains spring p/a cover data
Cleaned_Fall_Cover_Data_September_2019.xlsx - contains fall % cover data
Seed_Rain_Data_2019_FINAL.xlsx - contains seed rain data
Seed_Bank_2020_Cleaned.xlsx - contains seed bank data
SR_Traits_Dataset.xlsx - contains trait data (work in progress)
R Version: R 4.2.2
RStudio Version: 2023.03.0+386
Package Version:
Last updated: 16 June 2023
### --- Dataset importing and exporting -----
# Import in excel files
library(readxl)
# Export dataframes to excel files
library(writexl)
### ---- Data manipulation and visualization ----
# General data cleaning, manipulation, and visualization
library(tidyverse)
## ── Attaching core tidyverse packages ─────────────────── tidyverse 1.3.2.9000 ──
## ✔ dplyr 1.1.0 ✔ readr 2.1.3
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.0 ✔ tibble 3.1.8
## ✔ lubridate 1.9.1 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
#
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4
#
library(maditr)
##
## To drop variable use NULL: let(mtcars, am = NULL) %>% head()
##
##
## Attaching package: 'maditr'
##
## The following objects are masked from 'package:dplyr':
##
## between, coalesce, first, last
##
## The following object is masked from 'package:purrr':
##
## transpose
##
## The following object is masked from 'package:readr':
##
## cols
### ---- Model fitting and analysis ----
# Mixed-models
## Fits mixed models
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
## Functions to analyze significance and R^2 components of mixed-models
library(lmerTest)
##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## step
##
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
# Summary stats
##
library(mosaic)
## Registered S3 method overwritten by 'mosaic':
## method from
## fortify.SpatialPolygonsDataFrame ggplot2
##
## The 'mosaic' package masks several functions from core packages in order to add
## additional features. The original behavior of these functions should not be affected by this.
##
## Attaching package: 'mosaic'
##
## The following objects are masked from 'package:car':
##
## deltaMethod, logit
##
## The following object is masked from 'package:lmerTest':
##
## rand
##
## The following object is masked from 'package:lme4':
##
## factorize
##
## The following object is masked from 'package:Matrix':
##
## mean
##
## The following object is masked from 'package:permute':
##
## shuffle
##
## The following objects are masked from 'package:dplyr':
##
## count, do, tally
##
## The following object is masked from 'package:purrr':
##
## cross
##
## The following object is masked from 'package:ggplot2':
##
## stat
##
## The following objects are masked from 'package:stats':
##
## binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
## quantile, sd, t.test, var
##
## The following objects are masked from 'package:base':
##
## max, mean, min, prod, range, sample, sum
setwd("/Users/wynnekat/Desktop/Ph.D. - PFCA and TP Project Datasets")
spring_cover <- read_excel("/Users/wynnekat/Desktop/Ph.D. - PFCA and TP Project Datasets/Vegetative Cover Data 2019/Cleaned_Spring_Cover_Data_July_2019.xlsx")
fall_cover <- read_excel("/Users/wynnekat/Desktop/Ph.D. - PFCA and TP Project Datasets/Vegetative Cover Data 2019/Cleaned_Fall_Cover_Data_September_2019.xlsx")
seed_rain <- read_excel("/Users/wynnekat/Desktop/Ph.D. - PFCA and TP Project Datasets/Seed Rain 2019/Seed_Rain_Data_2019_FINAL.xlsx")
seed_bank <- read_excel("/Users/wynnekat/Desktop/Ph.D. - PFCA and TP Project Datasets/Seed Bank 2020/Seed_Bank_2020_Cleaned.xlsx")
traits <- read_excel("/Users/wynnekat/Desktop/Ph.D. - PFCA and TP Project Datasets/Traits and Species Info/SR_SB_VEG_Traits_Dataset.xlsx")
Study Sites
Our study sites were at Tucker Prairie Research Area (38○56’53.6” N, 91○59’40.0” W, Callaway County, MO) and at Prairie Fork Conservation Area (38○58’29.7” N, 91○44’03.3” W, Callaway County, MO). Tucker Prairie is a 59-hectare tract of unplowed North American tallgrass claypan prairie. Less than 0.5 % of intact tallgrass prairie ecosystem (i.e., never-been-plowed) remains in Missouri, and Tucker Prairie represents the last sizable remnant prairie in north central Missouri (Samson & Knopf, 1994). More than 250 species of plants inhabit Tucker Prairie, with representatives from 57 families and over 150 genera (Tropicos, n.d.). Native C4 grasses dominate the landscape including Sorghastrum nutans, Andropogon gerardii, Schizachyrium scoparium, and Sporobolus heterolepis (Drew, 1947; Rabinowitz & Rapp, 1980). Before 1957, Tucker Prairie was used for cattle grazing and occasional haying (Drew, 1947). From 1958 to 2002, Tucker Prairie was burned once every four years in the late winter or early spring (Rabinowitz & Rapp, 1980). Since 2002, Tucker Prairie has been managed on a 5-year burn rotation, where units are burned once in the late winter to early spring (Jan. – Mar.) and again 2-3 years later in the late summer to early fall (Aug. – Oct.).
Prairie Fork Conservation Area (PFCA) is over 450 hectares of former agricultural land being restored to tallgrass prairie and savanna ecosystems (Navarrete-Tindall et al., 2007; Newbold et al., 2019). From 2004 onwards, 16-25 hectares are newly seeded each year with native prairie species collected from Tucker Prairie and other nearby remnant prairies (Newbold et al., 2019). As a result, the plots within PFCA represent a chronosequence of reconstructed tallgrass prairies that are mainly comprised of Tucker Prairie descendants. Prior to seeding, glyphosate-resistant crops (i.e., corn and soybeans) are planted and harvested for at least three years to reduce the existing weed seed bank (Newbold et al., 2019). A seed mix containing an average of 179 species of native forbs, legumes, grasses, rushes, and sedges are broadcast seeded at a rate of 13.4 to 18.2 kg/ha once during the dormant season (Newbold et al., 2019). All mixes are comprised of at least 75% of the same species each year that are representative of nearby remnant prairies such as Tucker Prairie (Newbold et al., 2019). Reconstructions are managed using a 2-4-year burn schedule and annual spraying of invasive species (Newbold et al., 2019). For additional details on PFCA management, see Newbold et al. 2019. To capture changes in seed rain dynamics during the restoration process, we grouped our reconstructed sites into three categories, as defined by Newbold et al. (2019) as being compositionally different. We measured the seed rain in an old reconstruction (seeded in 2004), middle aged reconstruction (seeded in 2012 and 2013), and a young reconstruction (seeded in 2017).
Seed rain sampling
In May 2019, we deployed artificial turf grass seed traps (0.1 x 0.1 m) in Tucker Prairie and at each of the focal PFCA restorations. We used artificial turf grass traps instead of sticky traps due to their durability and resistance to freezing (Molau & Per Mølgaard, 1996; Bass, 2015). Turfgrass traps also do not lose effectiveness over time, like sticky traps, which are prone to contamination by soil, non-seed plant debris, and dead insects, making seed recovery difficult (Arruda et al., 2020). Furthermore, unlike funnel traps, they are easy to install and collect with minimal disturbance to surrounding vegetation and soil and do not fill with water (Schott, 1995). At each site, we randomly established ten, 5 m long transects, placed a minimum of 50 m apart from each other. We then placed traps at 1 m intervals along each transect and affixed them to the ground using ground staples (n = 5 traps per transect, 50 traps per site). We collected and replaced seed traps every 2-weeks during the growing season from May 31st to December 12th, 2019. Hereafter we refer to all captured diaspores as “seeds” despite catching a variety of fruit types (e.g., achene). After collection, all traps from a transect and sampling period were grouped together and sieved through a series of soil sieves (1 mm, 500 mm, 250 mm mesh). From these layers, we then picked out and identified seeds to the lowest taxonomic level possible using a reference collection of species inhabiting our study areas and identification guides (Coons et al., n.d.; A. C. Martin & Barkley, 1961).
Because of their extremely small size (0.3 – 0.5 mm; Yatskievych, 1999), we estimated the number of rush seeds (Juncus sp.) when there were over 200 seeds present in a sample. We first sieved the samples through 180 mm and 150 mm mesh soil sieves. Then we subsampled each sieved layer by calculating the average number of rush seeds per 1 cm2 area (n = 3) and multiplying that average by the total area covered by the sample layer. Lastly, we summed the number of estimated rush seeds per layer to calculate the total number of rush seeds in a sample.
Vegetation sampling
In 2019, at the same transects, we sampled species cover twice, once in the early summer (June - July) and again at peak biomass (September). During the first sampling period, we only recorded species presence. At peak biomass, we also sampled percent aerial species cover in a 1m2 area centered around each seed rain trap. We identified species according to Yatskievych (1999).
Seed bank sampling
In March 2020, before the growing season, we collected soil samples to assess the soil seed bank at each focal prairie. To compare the soil seed bank and the previous year’s seed rain, we collected 5 (10 x 10 x 10 cm3; 1000 cm3) soil cores approximately 0.5 m away from where we captured seed rain at each transect in 2019. We allowed the soil samples to dry before removing all non-seed plant material (e.g., roots and rhizomes) to ensure seedlings only germinated from seeds. We then homogenized all samples collected from a particular transect and subsampled ~1500 cm3 of soil to spread ~1 cm deep over sterile potting soil in plastic germination trays. Starting July 2020, we placed the trays (n = 40) in a greenhouse and watered them when dry. Amongst the 40 trays containing soil samples, we also placed an additional 10 trays containing only sterile potting soil to assess whether any contamination occurred from external sources. Once identifiable, we removed seedlings from trays. After germination ceased in July 2021, we placed all trays into a cold room for ~4 months to replicate conditions necessary to break dormancy for any still dormant seeds. Post-vernalization, we returned the trays to the greenhouse, stirred the soil samples, and monitored for any additional germination. We ended the study one year post-vernalization.
# Make the six letter species codes for the spring and fall cover uppercase
spring_cover[[3]] <- toupper(spring_cover[[3]])
fall_cover[[4]] <- toupper(fall_cover[[4]])
# Remove unnecessary columns from datasets
### Spring cover - Only keep Site, Transect, SPP6, Scientific Name
spring_cover <- spring_cover[,-c(5,6)]
### Fall cover - Keep everything
### Seed rain - Only keep Plot, Transect, Sampling Date, Week, SPP6, Number, Unknown
seed_rain <- seed_rain[, -c(7,8)]
### Seed bank - Keep everything
# Change names of columns to better reflect datasets (e.g. seed rain: number -> Number_Seeds)
colnames(seed_bank) <- c("Site", "Transect", "SPP6", "Number_Seedlings")
colnames(seed_rain) <- c("Site", "Transect", "Sampling_Date", "Week", "SPP6", "Number_Seeds", "Unknown")
# Change transect PFCA 2 - 10A to just 10 for the cover datasets
spring_cover <- spring_cover %>%
mutate_if(is.character,
str_replace_all, pattern = "10A", replacement = "10")
fall_cover <- fall_cover %>%
mutate_if(is.character,
str_replace_all, pattern = "10a", replacement = "10")
# Make Transect and Plot variables factors for all datasets
## Spring cover
spring_cover$Transect <- as.factor(spring_cover$Transect)
## Fall cover
fall_cover$Transect <- as.factor(fall_cover$Transect)
fall_cover$Plot <- as.factor(fall_cover$Plot)
## Seed Bank
seed_bank$Transect <- as.factor(seed_bank$Transect)
## Seed Rain
seed_rain$Transect <- as.factor(seed_rain$Transect)
####### Clean Spring
# 0.) Make a frequency column
spring_cover$Cover <- rep(1, nrow(spring_cover))
####### Clean Fall
# 0.) Make anything less than 1 = 1
fall_cover <- fall_cover %>%
mutate(Cover = ifelse(Cover < 1, 1, Cover))
# 1.) get a sum per species per transect of cover.
fall_sum <- fall_cover %>%
group_by(Site, Transect, SPP6, Scientific_Name) %>%
summarize(Cover = sum(Cover))
## `summarise()` has grouped output by 'Site', 'Transect', 'SPP6'. You can
## override using the `.groups` argument.
fall_sum
## # A tibble: 1,215 × 5
## # Groups: Site, Transect, SPP6 [1,214]
## Site Transect SPP6 Scientific_Name Cover
## <chr> <fct> <chr> <chr> <dbl>
## 1 PFCA 1 1 ACAVIR Acalypha virginica 11.5
## 2 PFCA 1 1 AMBART Ambrosia artemisiifolia 6
## 3 PFCA 1 1 BARVUL Barbarea vulgaris 3.5
## 4 PFCA 1 1 CHAFAS Chamaecrista fasciculata 388
## 5 PFCA 1 1 CONCAN Conyza canadensis 1
## 6 PFCA 1 1 ERISPP Erigeron spp. 2
## 7 PFCA 1 1 ERIVIL Eriochloa villosa 2
## 8 PFCA 1 1 EUPSER Eupatorium serotinum 3
## 9 PFCA 1 1 KUMSTI Kummerowia stipulacea 6
## 10 PFCA 1 1 LESCAP Lespedeza capitata 1
## # … with 1,205 more rows
####### Seed Rain
# 1.) get a sum per species per transect of cover.
seed_rain_sum <- seed_rain %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seeds = sum(Number_Seeds))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
# Clean in the same way you did the seed rain analysis
for(i in 1:nrow(seed_rain_sum)){
if(seed_rain_sum[i,3] == "AGAFAS"){seed_rain_sum[i,3] <- "AGASPP"}
if(seed_rain_sum[i,3] == "EUPSER"){seed_rain_sum[i,3] <- "EUPSPP"}
if(seed_rain_sum[i,3] == "GENPUB"){seed_rain_sum[i,3] <- "GENSPP"}
}
####### Seed Bank
# 1.) get a sum per species per transect of cover.
seed_bank_sum <- seed_bank %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seedlings = sum(Number_Seedlings))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
####### Remove Spring Unknowns
spring_cover <- filter(spring_cover, !grepl('UNK', SPP6))
# Looks good
list(unique(spring_cover$SPP6))
## [[1]]
## [1] "ACAVIR" "ACERUB" "ACHMIL" "AGRGIG" "AGRHYE" "ALOCAR" "AMBART" "AMBPSY"
## [9] "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB" "BAPAUS" "BAPBRA"
## [17] "BARVUL" "BIDARI" "BLECIL" "BROJAP" "CAMRAD" "CAPBUR" "CARSPP" "CARBIC"
## [25] "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN" "CORLAN" "CORPAL" "CORSPP"
## [33] "CORTRI" "CROSAG" "CYPACU" "CYPECH" "DALPUR" "DESILL" "DESSES" "DESSPP"
## [41] "DICLAN" "DIGISC" "ECHPAL" "ELETEN" "ELYVIR" "ERISPP" "ERYYUC" "EUPCOR"
## [49] "EUPSER" "EUTGYM" "FESARU" "FESPAR" "FRAVIR" "GALAPA" "GALOBT" "GENPUB"
## [57] "GERCAR" "HELFLE" "HELHEL" "HELMOL" "HELSPP" "HYPPUN" "JUNBRA" "JUNMAR"
## [65] "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "LACSER" "LEPVIR" "LESCAP" "LESCUN"
## [73] "LESVIR" "LIAPYC" "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP"
## [81] "MELSPP" "MONFIS" "MUHSPP" "MYOVER" "OENBIE" "OXADIL" "PARINT" "PENDIG"
## [89] "PLAVIR" "POAPRA" "POLSAN" "POLVER" "POTSIM" "PYCPIL" "PYCTEN" "RANABO"
## [97] "RATPIN" "RHUCOP" "ROSCAR" "ROSSPP" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI"
## [105] "SALAZU" "SCHSCO" "SENMAR" "SETFAB" "SETPAR" "SILANT" "SILINT" "SILLAC"
## [113] "SISSPP" "SOLALT" "SOLCAR" "SOLNEM" "SOLRIG" "SORNUT" "SPHOBT" "SPOHET"
## [121] "STRLEI" "SYMERI" "SYMLAN" "SYMLAT" "SYMNOV" "SYMOBL" "SYMOOL" "SYMORB"
## [129] "SYMPRA" "SYMSPP" "THLARV" "OENSPP" "TRAOHI" "TRIPER" "TRIREP" "TRISPP"
## [137] "ULMSPP" "SCISPP" "VERARV" "VERHEL" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT"
## [145] "ZIZAUR" "SCLTRI"
####### Remove Fall Unknowns
fall_sum <- filter(fall_sum, !grepl('UNK', SPP6))
# Looks good
list(unique(fall_sum$SPP6))
## [[1]]
## [1] "ACAVIR" "AMBART" "BARVUL" "CHAFAS" "CONCAN" "ERISPP" "ERIVIL" "EUPSER"
## [9] "KUMSTI" "LESCAP" "PENDIG" "PLAVIR" "RUDHIR" "SETFAB" "SOLALT" "STRLEI"
## [17] "SYMPIL" "ACHMIL" "CORTRI" "DALPUR" "DIGISC" "EUTGYM" "JUNSPP" "PYCTEN"
## [25] "RATPIN" "RUDSUB" "SCISPP" "SOLSPE" "SORNUT" "SYMERI" "SYMNOV" "SYMORB"
## [33] "TRICAM" "BROJAP" "HELMOL" "HYPPUN" "OENBIE" "SOLRIG" "SYMLAN" "AGRHYE"
## [41] "KUMSTR" "LESCUN" "LIAPYC" "PARINT" "SOLCAR" "TAROFF" "BAPALB" "ERASPE"
## [49] "JUNVIR" "LESVIR" "SYMPRA" "VERARV" "CERSPP" "ECHCRU" "OXADIL" "PLAMAJ"
## [57] "SILINT" "SILLAC" "VERSPP" "CORLAN" "CYPECH" "CYPSTR" "EUPCOR" "MEDLUP"
## [65] "SCHSCO" "BIDARI" "FESPAR" "MONFIS" "RANABO" "SOLNEM" "OENFIL" "ULMPUM"
## [73] "DESSPP" "ELYCAN" "AGATEN" "BAPBRA" "CARSPP" "CIRALT" "ERYYUC" "GENAND"
## [81] "KOEMAC" "LACSER" "PYCPIL" "RUMCRI" "SPOHET" "SYMOBL" "VERHEL" "BLECIL"
## [89] "BOLAST" "CYPACU" "SYMLAT" "SYMOOL" "ECHPAL" "MELSPP" "SYMANO" "SYMLAE"
## [97] "CORPAL" "ANDGER" "EUPALT" "POAPRA" "ALOCAR" "LIASPP" "LYTALA" "SPOCOM"
## [105] "SPILAC" "PANVIR" "SENMAR" "TRIFLA" "HELSPP" "ELYVIR" "HELHEL" "MYOVER"
## [113] "SALAZU" "TRAOHI" "BAPAUS" "LESPRO" "DESSES" "ROSCAR" "AMOCAN" "LEPDEN"
## [121] "CAMRAD" "RUBSPP" "ACERUB" "ZIZAUR" "DICLAN" "EUPPER" "GALOBT" "LYSLAN"
## [129] "POTSIM" "SETPAR" "SOLMIS" "VIOSAG" "AGAFAS" "ANTNEG" "CROSAG" "LINSUL"
## [137] "RHUGLA" "FRAVIR" "GENPUB" "MUHGLA" "JUNBRA" "ULMSPP" "POLVER" "TRISPP"
## [145] "CORSPP" "POLSAN" "SILANT" "TRIPER" "RHUCOP"
####### Remove Seed Rain Unknowns
# 1.) get a sum per species per transect of cover.
seed_rain_sum <- filter(seed_rain_sum, !grepl('UNK', SPP6))
seed_rain_sum <- filter(seed_rain_sum, !grepl('PANICUM?', SPP6))
seed_rain_sum <- filter(seed_rain_sum, !grepl('NONE', SPP6))
# Looks good
list(unique(seed_rain_sum$SPP6))
## [[1]]
## [1] "ACAVIR" "AMBART" "ANAMIN" "BOLAST" "CERSPP" "CHAFAS" "CHEALB"
## [8] "CONCAN" "DIGISC" "ERISPP" "ERIVIL" "EUPSPP" "KUMSPP" "LEPVIR"
## [15] "MEDLUP" "OXADIL" "PENDIG" "PLAVIR" "PYCTEN" "RUDHIR" "SETSPP"
## [22] "SOLALT" "SORNUT" "SYMSPP" "THLARV" "TRIPER" "VEROSP" "ALOCAR"
## [29] "CORTRI" "CYPECH" "DICLAN" "ERYYUC" "HYPPUN" "MONFIS" "MYOVER"
## [36] "OENBIE" "RATPIN" "SCIGEO" "SOLRIG" "SPHOBT" "STRLEI" "TRIFLA"
## [43] "ACHMIL" "BARVUL" "FESARU" "GERCAR" "LESCUN" "LIAPYC" "PANCAP"
## [50] "SYMNOV" "BIDARI" "CARBUS" "DESSPP" "ERASPE" "GALAPA" "LESVIR"
## [57] "POAPRA" "RUDSUB" "VERHAS" "VERSPP" "ANDGER" "CARFES" "ECHCRU"
## [64] "EREHIE" "PLAOCC" "SILINT" "AGASPP" "BROJAP" "CAPBUR" "CIRALT"
## [71] "CORLAN" "MELSPP" "MOLVER" "SCHSCO" "SCLTRI" "CYPSTR" "FESPAR"
## [78] "JUNSPP" "PYCPIL" "TAROFF" "AMATUB" "LESCAP" "OENFIL" "KOEMAC"
## [85] "DAUCAR" "LEUVUL" "VULOCT" "CYPACU" "HORPUS" "SCIPEN" "ECHPAL"
## [92] "SOLNEM" "CORPAL" "LINSUL" "DIAARM" "HELMOL" "PERLON" "RUMCRI"
## [99] "BLECIL" "LYTALA" "SPOHET" "ELESPP1" "SPOCOM" "CARCEP" "AGRHYE"
## [106] "GENSPP" "HELHEL" "TRAOHI" "BAPALB" "CARANN" "IPOHED" "AMOCAN"
## [113] "SALAZU" "CARBIC" "EUPCOR" "EUTGYM" "GALOBT" "RUBSPP" "SETPAR"
## [120] "ELESPP" "LYSLAN" "POLSPP" "HYPHIR" "LYCAME" "SILANT" "CARDSP"
## [127] "CROSAG" "LOBSPI" "VIOSAG"
####### Remove Seed Bank Unknowns
seed_bank_sum<- filter(seed_bank_sum, !grepl('UNK', SPP6))
# Looks good
list(sort(unique(seed_bank_sum$SPP6)))
## [[1]]
## [1] "ABUTHE" "ACAVIR" "ACHMIL" "AGRHYE" "ALOCAR" "AMATUB" "AMBART" "ANDGER"
## [9] "BARVUL" "CARFES" "CARHIR" "CARPAR" "CERGLO" "CHAFAS" "CIRALT" "CONCAN"
## [17] "CORLAN" "CROSAG" "CYPACU" "CYPSPP" "CYPSTR" "DESSPP" "DICLAN" "DIGISC"
## [25] "DIGSAN" "DIPLAC" "ERASPE" "ERIANN" "ERISTR" "EUPHUM" "EUPSER" "EUTGYM"
## [33] "GERCAR" "HORPUS" "HYPPUN" "IPOHED" "JUNINT" "JUNMAR" "KUMSTI" "KUMSTR"
## [41] "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LYSLAN" "MEDLUP" "MELSPP" "MOLVER"
## [49] "OENSPP" "OXADIL" "PANCAP" "PANDIC" "PENDIG" "PERLON" "PLAOCC" "PLAPUS"
## [57] "PLAVIR" "POACHA" "POAPRA" "POPDEL" "POROLE" "PYCPIL" "PYCTEN" "RATPIN"
## [65] "RUDHIR" "RUMCRI" "SCHSCO" "SETFAB" "SETGLA" "SILANT" "SOLALT" "SOLCAR"
## [73] "SORNUT" "SPHOBT" "STRLEI" "SYMLAE" "SYMPIL" "TAROFF" "THLARV" "TOXRAD"
## [81] "TRAOHI" "TRIFLA" "TRIPER" "TRIREP" "VERHAS" "VERPER" "VERTHA"
# Full join the two cover data sets
cover_merged <- full_join(spring_cover, fall_sum, by = c("Site", "Transect", "SPP6", "Scientific_Name", "Cover")) %>%
mutate_all(~replace(., is.na(.), 0))
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
# Make sure there is no weirdness in the join
#View(cover_merged)
### Group by and then select the SPP6 with the max cover between the two datasets
### Since Spring only had presence/absence we don't want to inflate the fall cover artificially if we saw that species again
cover_merged_max <- cover_merged %>%
group_by(Site, Transect, SPP6, Scientific_Name) %>%
summarise(Cover=max(Cover))
## `summarise()` has grouped output by 'Site', 'Transect', 'SPP6'. You can
## override using the `.groups` argument.
## Not excatly sure what to do about SYMSPP and KUMSPP (things we couldn't ID in the spring but could in the fall)
### Going to drop SYMSPP from the dataset since every plot that had SYMSPP in the spring had identified Symphyotrichum in the fall
cover_merged_max <- filter(cover_merged_max , !grepl('SYMSPP', SPP6))
### Checked KUMSPP only PFCA 2 - 3 and PFCA 2 - 9 did not see Kummerowia again in the fall; will drop KUMSPP for all other plots
### Removes all the PFCA 1 KUMSPPs
cover_merged_max <- cover_merged_max[!(cover_merged_max$SPP6%in%"KUMSPP" & cover_merged_max$Site%in% "PFCA 1"),]
### Removes all PFCA 2 KUMSPPs that were not in PFCA 2 - 3 and PFCA 2- 9
cover_merged_max <- cover_merged_max[!(cover_merged_max$SPP6%in%"KUMSPP" & cover_merged_max$Site%in% "PFCA 2" & cover_merged_max$Transect%in% c(2,4,5,8)),]
### Check all the species codes
sort(unique(cover_merged_max$SPP6))
## [1] "ACAVIR" "ACERUB" "ACHMIL" "AGAFAS" "AGATEN" "AGRGIG" "AGRHYE" "ALOCAR"
## [9] "AMBART" "AMBPSY" "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB"
## [17] "BAPAUS" "BAPBRA" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP" "CAMRAD"
## [25] "CAPBUR" "CARBIC" "CARSPP" "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN"
## [33] "CORLAN" "CORPAL" "CORSPP" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR"
## [41] "DALPUR" "DESILL" "DESSES" "DESSPP" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL"
## [49] "ELETEN" "ELYCAN" "ELYVIR" "ERASPE" "ERISPP" "ERIVIL" "ERYYUC" "EUPALT"
## [57] "EUPCOR" "EUPPER" "EUPSER" "EUTGYM" "FESARU" "FESPAR" "FRAVIR" "GALAPA"
## [65] "GALOBT" "GENAND" "GENPUB" "GERCAR" "HELFLE" "HELHEL" "HELMOL" "HELSPP"
## [73] "HYPPUN" "JUNBRA" "JUNMAR" "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "KUMSTI"
## [81] "KUMSTR" "LACSER" "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LESPRO" "LESVIR"
## [89] "LIAPYC" "LIASPP" "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP"
## [97] "MELSPP" "MONFIS" "MUHGLA" "MUHSPP" "MYOVER" "OENBIE" "OENFIL" "OENSPP"
## [105] "OXADIL" "PANVIR" "PARINT" "PENDIG" "PLAMAJ" "PLAVIR" "POAPRA" "POLSAN"
## [113] "POLVER" "POTSIM" "PYCPIL" "PYCTEN" "RANABO" "RATPIN" "RHUCOP" "RHUGLA"
## [121] "ROSCAR" "ROSSPP" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU" "SCHSCO"
## [129] "SCISPP" "SCLTRI" "SENMAR" "SETFAB" "SETPAR" "SILANT" "SILINT" "SILLAC"
## [137] "SISSPP" "SOLALT" "SOLCAR" "SOLMIS" "SOLNEM" "SOLRIG" "SOLSPE" "SORNUT"
## [145] "SPHOBT" "SPILAC" "SPOCOM" "SPOHET" "STRLEI" "SYMANO" "SYMERI" "SYMLAE"
## [153] "SYMLAN" "SYMLAT" "SYMNOV" "SYMOBL" "SYMOOL" "SYMORB" "SYMPIL" "SYMPRA"
## [161] "TAROFF" "THLARV" "TRAOHI" "TRICAM" "TRIFLA" "TRIPER" "TRIREP" "TRISPP"
## [169] "ULMPUM" "ULMSPP" "VERARV" "VERHEL" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT"
## [177] "ZIZAUR"
### Cover grouping so that it works with the seed rain dataset
# Code to lump certain species together
cover_to_merge_rain <- cover_merged_max
for(i in 1:nrow(cover_to_merge_rain)){
# Group all Symphyotrichum except SYMNOV
## SYMPIL
## SYMLAN
## SYMLAT
## SYMERI
## SYMANO
## SYMLAE
## SYMPRA
## SYMOBL
## SYMOOL
if(cover_to_merge_rain[i,3] == "SYMPIL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMLAN"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMLAT"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMERI"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMANO"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMLAE"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMPRA"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMOBL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMOOL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
# Group all Carex
## CARBIC
if(cover_to_merge_rain[i,3] == "CARBIC"){cover_to_merge_rain[i,3] <- "CARSPP"}
# Group all Kummerowia
## KUMSTI
## KUMSTR
if(cover_to_merge_rain[i,3] == "KUMSTI"){cover_to_merge_rain[i,3] <- "KUMSPP"}
if(cover_to_merge_rain[i,3] == "KUMSTR"){cover_to_merge_rain[i,3] <- "KUMSPP"}
# Group all Juncus
## JUNMAR
## JUNBBRA
if(cover_to_merge_rain[i,3] == "JUNMAR"){cover_to_merge_rain[i,3] <- "JUNSPP"}
if(cover_to_merge_rain[i,3] == "JUNBRA"){cover_to_merge_rain[i,3] <- "JUNSPP"}
# Group all Agalinis
## AGATEN
## AGAFAS
if(cover_to_merge_rain[i,3] == "AGATEN"){cover_to_merge_rain[i,3] <- "AGASPP"}
if(cover_to_merge_rain[i,3] == "AGAFAS"){cover_to_merge_rain[i,3] <- "AGASPP"}
# Group all Eleocharis
## ELETEN
if(cover_to_merge_rain[i,3] == "ELETEN"){cover_to_merge_rain[i,3] <- "ELESPP"}
# Group all Gentian
## GENPUB
## GENAND
if(cover_to_merge_rain[i,3] == "GENPUB"){cover_to_merge_rain[i,3] <- "GENSPP"}
if(cover_to_merge_rain[i,3] == "GENAND"){cover_to_merge_rain[i,3] <- "GENSPP"}
# Group all Setaria except Setaria parviflora
## SETGLA
## SETFAB
if(cover_to_merge_rain[i,3] == "SETGLA"){cover_to_merge_rain[i,3] <- "SETSPP"}
if(cover_to_merge_rain[i,3] == "SETFAB"){cover_to_merge_rain[i,3] <- "SETSPP"}
# Group all Veronica
## VERARV
if(cover_to_merge_rain[i,3] == "VERARV"){cover_to_merge_rain[i,3] <- "VEROSP"}
# Group all Desmodium
## DESSES
if(cover_to_merge_rain[i,3] == "DESSES"){cover_to_merge_rain[i,3] <- "DESSPP"}
# Group all Eupatorium
## EUPALT
## EUPSER
## EUPPER
if(cover_to_merge_rain[i,3] == "EUPALT"){cover_to_merge_rain[i,3] <- "EUPSPP"}
if(cover_to_merge_rain[i,3] == "EUPSER"){cover_to_merge_rain[i,3] <- "EUPSPP"}
if(cover_to_merge_rain[i,3] == "EUPPER"){cover_to_merge_rain[i,3] <- "EUPSPP"}
# Group all Polygala
## POLVER
## POLSAN
if(cover_to_merge_rain[i,3] == "POLVER"){cover_to_merge_rain[i,3] <- "POLSPP"}
if(cover_to_merge_rain[i,3] == "POLSAN"){cover_to_merge_rain[i,3] <- "POLSPP"}
}
### Check that it changed the codes
sort(unique(cover_to_merge_rain$SPP6))
## [1] "ACAVIR" "ACERUB" "ACHMIL" "AGASPP" "AGRGIG" "AGRHYE" "ALOCAR" "AMBART"
## [9] "AMBPSY" "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB" "BAPAUS"
## [17] "BAPBRA" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP" "CAMRAD" "CAPBUR"
## [25] "CARSPP" "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN" "CORLAN" "CORPAL"
## [33] "CORSPP" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR" "DALPUR" "DESILL"
## [41] "DESSPP" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL" "ELESPP" "ELYCAN" "ELYVIR"
## [49] "ERASPE" "ERISPP" "ERIVIL" "ERYYUC" "EUPCOR" "EUPSPP" "EUTGYM" "FESARU"
## [57] "FESPAR" "FRAVIR" "GALAPA" "GALOBT" "GENSPP" "GERCAR" "HELFLE" "HELHEL"
## [65] "HELMOL" "HELSPP" "HYPPUN" "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "LACSER"
## [73] "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LESPRO" "LESVIR" "LIAPYC" "LIASPP"
## [81] "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP" "MELSPP" "MONFIS"
## [89] "MUHGLA" "MUHSPP" "MYOVER" "OENBIE" "OENFIL" "OENSPP" "OXADIL" "PANVIR"
## [97] "PARINT" "PENDIG" "PLAMAJ" "PLAVIR" "POAPRA" "POLSPP" "POTSIM" "PYCPIL"
## [105] "PYCTEN" "RANABO" "RATPIN" "RHUCOP" "RHUGLA" "ROSCAR" "ROSSPP" "RUBSPP"
## [113] "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU" "SCHSCO" "SCISPP" "SCLTRI" "SENMAR"
## [121] "SETPAR" "SETSPP" "SILANT" "SILINT" "SILLAC" "SISSPP" "SOLALT" "SOLCAR"
## [129] "SOLMIS" "SOLNEM" "SOLRIG" "SOLSPE" "SORNUT" "SPHOBT" "SPILAC" "SPOCOM"
## [137] "SPOHET" "STRLEI" "SYMNOV" "SYMORB" "SYMSPP" "TAROFF" "THLARV" "TRAOHI"
## [145] "TRICAM" "TRIFLA" "TRIPER" "TRIREP" "TRISPP" "ULMPUM" "ULMSPP" "VERHEL"
## [153] "VEROSP" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT" "ZIZAUR"
### Sum everything by species by transect again
cover_to_merge_rain <- cover_to_merge_rain %>%
group_by(Site, Transect, SPP6, Scientific_Name) %>%
summarize(Cover = sum(Cover))
## `summarise()` has grouped output by 'Site', 'Transect', 'SPP6'. You can
## override using the `.groups` argument.
### Seed grouping so that it works with the cover dataset
seed_rain_names <- left_join(seed_rain_sum, traits)
## Joining with `by = join_by(SPP6)`
rain_to_merge_cover <- seed_rain_names[,c(1,2,3,4,5)]
for(i in 1:nrow(rain_to_merge_cover)){
# Group all Carex
## CARBUS
## CARBIC
## CARANN
## CARFES
## CARCEP
if(rain_to_merge_cover[i,3] == "CARBUS"){rain_to_merge_cover[i,3] <- "CARSPP"}
if(rain_to_merge_cover[i,3] == "CARBIC"){rain_to_merge_cover[i,3] <- "CARSPP"}
if(rain_to_merge_cover[i,3] == "CARANN"){rain_to_merge_cover[i,3] <- "CARSPP"}
if(rain_to_merge_cover[i,3] == "CARFES"){rain_to_merge_cover[i,3] <- "CARSPP"}
if(rain_to_merge_cover[i,3] == "CARCEP"){rain_to_merge_cover[i,3] <- "CARSPP"}
# Group all Eleocharis
## ELESPP1
if(rain_to_merge_cover[i,3] == "ELESPP1"){rain_to_merge_cover[i,3] <- "ELESPP"}
# Group all Scirpus
## ELESPP1
if(rain_to_merge_cover[i,3] == "SCIGEO"){rain_to_merge_cover[i,3] <- "SCISPP"}
if(rain_to_merge_cover[i,3] == "SCIPEN"){rain_to_merge_cover[i,3] <- "SCISPP"}
}
### Sum everything by species by transect again
rain_to_merge_cover <- rain_to_merge_cover %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seeds = sum(Number_Seeds))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
# Looks good
sort(unique(rain_to_merge_cover$SPP6))
## [1] "ACAVIR" "ACHMIL" "AGASPP" "AGRHYE" "ALOCAR" "AMATUB" "AMBART" "AMOCAN"
## [9] "ANAMIN" "ANDGER" "BAPALB" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP"
## [17] "CAPBUR" "CARDSP" "CARSPP" "CERSPP" "CHAFAS" "CHEALB" "CIRALT" "CONCAN"
## [25] "CORLAN" "CORPAL" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR" "DAUCAR"
## [33] "DESSPP" "DIAARM" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL" "ELESPP" "ERASPE"
## [41] "EREHIE" "ERISPP" "ERIVIL" "ERYYUC" "EUPCOR" "EUPSPP" "EUTGYM" "FESARU"
## [49] "FESPAR" "GALAPA" "GALOBT" "GENSPP" "GERCAR" "HELHEL" "HELMOL" "HORPUS"
## [57] "HYPHIR" "HYPPUN" "IPOHED" "JUNSPP" "KOEMAC" "KUMSPP" "LEPVIR" "LESCAP"
## [65] "LESCUN" "LESVIR" "LEUVUL" "LIAPYC" "LINSUL" "LOBSPI" "LYCAME" "LYSLAN"
## [73] "LYTALA" "MEDLUP" "MELSPP" "MOLVER" "MONFIS" "MYOVER" "OENBIE" "OENFIL"
## [81] "OXADIL" "PANCAP" "PENDIG" "PERLON" "PLAOCC" "PLAVIR" "POAPRA" "POLSPP"
## [89] "PYCPIL" "PYCTEN" "RATPIN" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU"
## [97] "SCHSCO" "SCISPP" "SCLTRI" "SETPAR" "SETSPP" "SILANT" "SILINT" "SOLALT"
## [105] "SOLNEM" "SOLRIG" "SORNUT" "SPHOBT" "SPOCOM" "SPOHET" "STRLEI" "SYMNOV"
## [113] "SYMSPP" "TAROFF" "THLARV" "TRAOHI" "TRIFLA" "TRIPER" "VERHAS" "VEROSP"
## [121] "VERSPP" "VIOSAG" "VULOCT"
# Drop the scientific names they are causing trouble with the join
cover_to_merge_rain <- cover_to_merge_rain[, -4]
cover_rain <- full_join(cover_to_merge_rain, rain_to_merge_cover) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:
# Drop TP transect 7 because of missing data
cover_rain <- cover_rain[!(cover_rain$Site%in%"TP" & cover_rain$Transect%in% "7"),]
### Cover grouping so that it works with the seed bank dataset
# Code to lump certain species together
cover_to_merge_bank <- cover_merged_max
for(i in 1:nrow(cover_to_merge_bank)){
# Group all Desmodium together
## DESSES
if(cover_to_merge_bank[i,3] == "DESSES"){cover_to_merge_bank[i,3] <- "DESSPP"}
# Group all Cyperus together
## CYPSTR
## CYPECH
if(cover_to_merge_bank[i,3] == "CYPSTR"){cover_to_merge_bank[i,3] <- "CYPSPP"}
if(cover_to_merge_bank[i,3] == "CYPECH"){cover_to_merge_bank[i,3] <- "CYPSPP"}
# Group all Oenothera together
## OENFIL
## OENBIE
if(cover_to_merge_bank[i,3] == "OENFIL"){cover_to_merge_bank[i,3] <- "OENSPP"}
if(cover_to_merge_bank[i,3] == "OENBIE"){cover_to_merge_bank[i,3] <- "OENSPP"}
}
cover_to_merge_bank <- cover_to_merge_bank %>%
group_by(Site, Transect, SPP6) %>%
summarize(Cover = sum(Cover))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
### Seed Bank grouping so that it works with the cover dataset
bank_to_merge_cover <- seed_bank_sum
# Code to lump certain species together
for(i in 1:nrow(bank_to_merge_cover)){
# Group all Cerastium together
## CERGLO
if(bank_to_merge_cover[i,3] == "CERGLO"){bank_to_merge_cover[i,3] <- "CERSPP"}
# Group all Cyperus together
## DESSES
if(bank_to_merge_cover[i,3] == "CYPSTR"){bank_to_merge_cover[i,3] <- "CYPSPP"}
# Group all the Carex together
if(bank_to_merge_cover[i,3] == "CARFES"){bank_to_merge_cover[i,3] <- "CARSPP"}
# Group all Erigeron together
## ERIANN
## ERISTR
if(bank_to_merge_cover[i,3] == "ERIANN"){bank_to_merge_cover[i,3] <- "ERISPP"}
if(bank_to_merge_cover[i,3] == "ERISTR"){bank_to_merge_cover[i,3] <- "ERISPP"}
# Change Juncus interior to Juncus sp. (since we were able to ID JUNMAR in the cover but not JUNINT)
if(bank_to_merge_cover[i,3] == "JUNINT"){bank_to_merge_cover[i,3] <- "JUNSPP"}
}
bank_to_merge_cover <- bank_to_merge_cover %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seedlings = sum(Number_Seedlings))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
# Drop the scientific names they are causing trouble with the join
cover_bank <- full_join(cover_to_merge_bank, bank_to_merge_cover) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:
# Drop TP transect 7 because of missing fall cover data
cover_bank <- cover_bank[!(cover_bank$Site%in%"TP" & cover_bank$Transect%in% "7"),]
### Seed Bank grouping so that it works with the cover dataset
bank_to_merge_rain <- seed_bank_sum
# Code to lump certain species together
for(i in 1:nrow(bank_to_merge_rain)){
# Group all Cerastium together
## CERGLO
if(bank_to_merge_rain[i,3] == "CERGLO"){bank_to_merge_rain[i,3] <- "CERSPP"}
# Group all Cyperus together
## CYPSTR
if(bank_to_merge_rain[i,3] == "CYPSTR"){bank_to_merge_rain[i,3] <- "CYPSPP"}
# Group all Erigeron together
## ERIANN
## ERISTR
if(bank_to_merge_rain[i,3] == "ERIANN"){bank_to_merge_rain[i,3] <- "ERISPP"}
if(bank_to_merge_rain[i,3] == "ERISTR"){bank_to_merge_rain[i,3] <- "ERISPP"}
# Group all the Juncus together
## JUNINT
## JUNMAR
if(bank_to_merge_rain[i,3] == "JUNINT"){bank_to_merge_rain[i,3] <- "JUNSPP"}
if(bank_to_merge_rain[i,3] == "JUNMAR"){bank_to_merge_rain[i,3] <- "JUNSPP"}
# Group all the Setaria together
if(bank_to_merge_rain[i,3] == "SETFAB"){bank_to_merge_rain[i,3] <- "SETSPP"}
if(bank_to_merge_rain[i,3] == "SETGLA"){bank_to_merge_rain[i,3] <- "SETSPP"}
# Group all the Symphyotrichum together
## SYMPIL
## SYMLAE
if(bank_to_merge_rain[i,3] == "SYMPIL"){bank_to_merge_rain[i,3] <- "SYMSPP"}
if(bank_to_merge_rain[i,3] == "SYMLAE"){bank_to_merge_rain[i,3] <- "SYMSPP"}
if(bank_to_merge_rain[i,3] == "SIMPIL"){bank_to_merge_rain[i,3] <- "SYMSPP"}
# Group all the Kummerowia together
## KUMSTR
## KUMSTI
if(bank_to_merge_rain[i,3] == "KUMSTR"){bank_to_merge_rain[i,3] <- "KUMSPP"}
if(bank_to_merge_rain[i,3] == "KUMSTI"){bank_to_merge_rain[i,3] <- "KUMSPP"}
# Group all the Veronica together
## VERPER
if(bank_to_merge_rain[i,3] == "VERPER"){bank_to_merge_rain[i,3] <- "VEROSP"}
# Group all the Eupatorium together
## EUPSER
if(bank_to_merge_rain[i,3] == "EUPSER"){bank_to_merge_rain[i,3] <- "EUPSPP"}
# Merge all the Cardamine together
## CARHIR
## CARPAR
if(bank_to_merge_rain[i,3] == "CARHIR"){bank_to_merge_rain[i,3] <- "CARDSP"}
if(bank_to_merge_rain[i,3] == "CARPAR"){bank_to_merge_rain[i,3] <- "CARDSP"}
}
bank_to_merge_rain <- bank_to_merge_rain %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seedlings = sum(Number_Seedlings))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
rain_to_merge_bank <- seed_rain_sum
for(i in 1:nrow(rain_to_merge_bank)){
# Group all Cyperus together
## CYPSTR
## CYPECH
if(rain_to_merge_bank[i,3] == "CYPSTR"){rain_to_merge_bank[i,3] <- "CYPSPP"}
if(rain_to_merge_bank[i,3] == "CYPECH"){rain_to_merge_bank[i,3] <- "CYPSPP"}
# Group all Oenothera together
## OENFIL
## OENBIE
if(rain_to_merge_bank[i,3] == "OENFIL"){rain_to_merge_bank[i,3] <- "OENSPP"}
if(rain_to_merge_bank[i,3] == "OENBIE"){rain_to_merge_bank[i,3] <- "OENSPP"}
}
rain_to_merge_bank <- rain_to_merge_bank %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seeds = sum(Number_Seeds))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
rain_bank <- full_join(bank_to_merge_rain, rain_to_merge_bank) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:
cover_rain_sum <- cover_rain %>%
group_by(SPP6) %>%
summarize(tot.cover = sum(Cover), tot.seeds = sum(Number_Seeds))
Wasn’t in the cover at all but we caught seeds
Didn’t catch seeds but was in the cover
Notable (> 10% cover observed)
Minor ( 1% < cover < 10% observed)
Trace (cover 1% observed)
## # A tibble: 2,750 × 5
## # Groups: Site, Transect, SPP6 [1,809]
## Site Transect SPP6 Abundance Type
## <chr> <fct> <chr> <dbl> <chr>
## 1 PFCA 1 1 ACAVIR 11.5 Aboveground
## 2 PFCA 1 1 ALOCAR 1 Aboveground
## 3 PFCA 1 1 AMBART 6 Aboveground
## 4 PFCA 1 1 BARVUL 3.5 Aboveground
## 5 PFCA 1 1 BROJAP 1 Aboveground
## 6 PFCA 1 1 CERSPP 1 Aboveground
## 7 PFCA 1 1 CHAFAS 388 Aboveground
## 8 PFCA 1 1 CONCAN 1 Aboveground
## 9 PFCA 1 1 DESILL 1 Aboveground
## 10 PFCA 1 1 ERISPP 2 Aboveground
## # … with 2,740 more rows
# Code obtained and modified from Lauren and Anu's seed rain vs. vegetation paper
##### run dissimilarity loops
Cover_rain_dis_results <-matrix(nrow=0,ncol=7)
sites <-as.vector(unique(cover_rain_dis$Site))
for(s in 1:length(sites)){
temp <-subset(cover_rain_dis, Site==sites[s])
temp$SPP6 <- factor(as.character(temp$SPP6))
transects <-as.vector(unique(temp$Transect))
#for(b in 1:length(blocks)){
#temp_b<-subset(temp, block==blocks[b])
#plots<-as.vector(unique(temp_b$plot))
for(t in 1:length(transects)){
temp_t <- subset(temp, Transect==transects[t])
#wide form - with a row for above and belowground
plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
#relative abundances
trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"tot"))
names(trans_tot)[1]<-"trt"
trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
names(trans_pa)[1]<-"trt"
trans_int <- cbind(plot_cast$Type,trunc(plot_cast[,-1])) #integers only
names(trans_int)[1]<-"trt"
trans_int_norm <- cbind(plot_cast$Type,trunc(decostand(plot_cast[,-1], "normalize")*100)) #integers only and normalize
names(trans_int)[1]<-"trt"
distance_bray<-vegdist(trans_tot[,-1], method = "bray")
distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")
distance_kulczynski<-vegdist(trans_tot[,-1], method = "kulczynski")
distance_cao<-vegdist(trans_int_norm[,-1], method = "cao")
distance_morisita<-vegdist(trans_int[,-1], method = "morisita")
new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1],
distance_cao[1], distance_morisita[1], distance_kulczynski[1])
Cover_rain_dis_results <-rbind(Cover_rain_dis_results, new.row)
}
print(s/length(sites))
}
## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1
row.names(Cover_rain_dis_results)<-NULL
Cover_rain_dis_results<- data.frame(Cover_rain_dis_results)
names(Cover_rain_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard",
"dissimilarity_cao", "dissimilarity_morisita", "dissimilarity_kulczynski")
##
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = Cover_rain_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.267522 -0.063391 -0.002319 0.088438 0.190470
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.46881 0.03621 12.949 6.58e-15 ***
## SitePFCA 2 0.16217 0.05120 3.167 0.003186 **
## SitePFCA 3 0.21871 0.05120 4.272 0.000141 ***
## SiteTP 0.21184 0.05261 4.027 0.000289 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1145 on 35 degrees of freedom
## Multiple R-squared: 0.4014, Adjusted R-squared: 0.3501
## F-statistic: 7.823 on 3 and 35 DF, p-value: 0.0003994
## Anova Table (Type III tests)
##
## Response: dissimilarity_bray
## Sum Sq Df F value Pr(>F)
## (Intercept) 2.19783 1 167.6676 6.577e-15 ***
## Site 0.30764 3 7.8231 0.0003994 ***
## Residuals 0.45879 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 95% confidence intervals
emmeans::emmeans(Cover_rain.mod.bray, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## PFCA 1 0.469 0.0362 35 0.395 0.542
## PFCA 2 0.631 0.0362 35 0.557 0.704
## PFCA 3 0.688 0.0362 35 0.614 0.761
## TP 0.681 0.0382 35 0.603 0.758
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.472, 0.628, 0.687, 0.681)
lower <- c(0.396, 0.552, 0.611, 0.601)
upper <- c(0.548, 0.704, 0.763, 0.761)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Cover_rain_dissimilarity <- data.frame(site, response, lower, upper)
## 95% confidence intervals
emmeans::emmeans(Cover_rain.mod.jac, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## PFCA 1 0.480 0.0193 35 0.440 0.519
## PFCA 2 0.462 0.0193 35 0.423 0.502
## PFCA 3 0.589 0.0193 35 0.550 0.628
## TP 0.542 0.0204 35 0.500 0.583
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.480, 0.462, 0.589, 0.542)
lower <- c(0.440, .423, 0.550, 0.500)
upper <- c(0.519, 0.502, 0.628, 0.583)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Cover_rain_jac_dissimilarity <- data.frame(site, response, lower, upper)
cover_rain_dis_wide <- cover_rain_dis %>%
spread(SPP6, Abundance) %>%
mutate_all(~replace(., is.na(.), 0))
## `mutate_all()` ignored the following grouping variables:
## • Columns `Site`, `Transect`
## ℹ Use `mutate_at(df, vars(-group_cols()), myoperation)` to silence the message.
cover_rain_dis_long <- cover_rain_dis_wide %>%
gather(key = "SPP6", value = "Abundance", 4:173) %>%
spread(key = Type, Abundance)
# Get rid of columns that both have zeros
cover_rain_dis_long_unique <- cover_rain_dis_long[!(cover_rain_dis_long$Aboveground == 0 & cover_rain_dis_long$Seed_Rain == 0),]
# 0 = no seeds/seedlings and 1 = seeds/seedlings
cover_rain_dis_long_unique <- cover_rain_dis_long_unique %>%
mutate(Aboveground = if_else(Aboveground > 0, 1, 0)) %>%
mutate(Seed_Rain = if_else(Seed_Rain > 0, 1, 0))
# Make a new column that adds both the seed rain and aboveground presence/absence columns. Then remove rows that = 2 (shares both species)
cover_rain_dis_long_unique$shared <- cover_rain_dis_long_unique$Aboveground + cover_rain_dis_long_unique$Seed_Rain
cover_rain_dis_long_unique2 <- cover_rain_dis_long_unique[!(cover_rain_dis_long_unique$shared == 2),]
cover_rain_dis_long_unique2 <- cover_rain_dis_long_unique2[, -6]
# Count unique taxa per site and transect for the seed bank compared to the seed rain
Unique_species_aboveground_com_rain <- cover_rain_dis_long_unique2[, -5]
Unique_species_aboveground_com_rain <- Unique_species_aboveground_com_rain %>%
filter(Aboveground != 0 ) %>%
group_by(Site, Transect) %>%
summarize(New_Aboveground = length(unique(SPP6)))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
# Count unique taxa per site and transect for the seed rain compared to the seed bank
Unique_species_rain_com_aboveground <- cover_rain_dis_long_unique2[, -4]
Unique_species_rain_com_aboveground <- Unique_species_rain_com_aboveground %>%
filter(Seed_Rain != 0 ) %>%
group_by(Site, Transect) %>%
summarize(New_Rain = length(unique(SPP6)))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
# Count the number of shared taxa between communities
cover_rain_dis_long_shared <- cover_rain_dis_long_unique %>%
filter(shared == 2) %>%
group_by(Site, Transect) %>%
summarize(Shared = length(unique(SPP6)))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
## Merge them together
Species_counts_cover_rain <- full_join(Unique_species_aboveground_com_rain, Unique_species_rain_com_aboveground) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:
## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 3 remaining warnings.
Species_counts_cover_rain <- full_join(Species_counts_cover_rain, cover_rain_dis_long_shared) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:
## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 3 remaining warnings.
Species_counts_cover_rain$Tot_Richness <- Species_counts_cover_rain$New_Aboveground + Species_counts_cover_rain$New_Rain + Species_counts_cover_rain$Shared
# Count the number of species in the seed rain for each transect
Species_counts_cover_rain$Tot_Rain <- Species_counts_cover_rain$New_Rain + Species_counts_cover_rain$Shared
# Count the number of species in the seed bank for each transect
Species_counts_cover_rain$Tot_Aboveground <- Species_counts_cover_rain$New_Aboveground + Species_counts_cover_rain$Shared
## Immigrants
Species_counts_cover_rain$Prop_Immigrant <- Species_counts_cover_rain$New_Rain / Species_counts_cover_rain$Tot_Aboveground
favstats(Species_counts_cover_rain$Prop_Immigrant~Species_counts_cover_rain$Site)
## Species_counts_cover_rain$Site min Q1 median Q3
## 1 PFCA 1 0.2325581 0.2454268 0.3009259 0.4000000
## 2 PFCA 2 0.1621622 0.2231378 0.2889714 0.3224343
## 3 PFCA 3 0.1081081 0.2181818 0.3042017 0.3495879
## 4 TP 0.1851852 0.2162162 0.2592593 0.3103448
## max mean sd n missing
## 1 0.6190476 0.3420640 0.12881067 10 0
## 2 0.3555556 0.2738325 0.06337656 10 0
## 3 0.4736842 0.2832612 0.11509477 10 0
## 4 0.5161290 0.2849422 0.10421660 9 0
# No difference in proportion of immigrant species
## *Note* better modeled as a binomial regression?
Species_counts_cover_rain.mod <- lm(Prop_Immigrant~Site, data =Species_counts_cover_rain)
summary(Species_counts_cover_rain.mod)
##
## Call:
## lm(formula = Prop_Immigrant ~ Site, data = Species_counts_cover_rain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.17515 -0.07599 -0.01296 0.05072 0.27698
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.34206 0.03345 10.226 4.7e-12 ***
## SitePFCA 2 -0.06823 0.04730 -1.442 0.158
## SitePFCA 3 -0.05880 0.04730 -1.243 0.222
## SiteTP -0.05712 0.04860 -1.175 0.248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1058 on 35 degrees of freedom
## Multiple R-squared: 0.06862, Adjusted R-squared: -0.01121
## F-statistic: 0.8596 on 3 and 35 DF, p-value: 0.4711
Anova(Species_counts_cover_rain.mod, type = "III")
## Anova Table (Type III tests)
##
## Response: Prop_Immigrant
## Sum Sq Df F value Pr(>F)
## (Intercept) 1.17008 1 104.5809 4.703e-12 ***
## Site 0.02885 3 0.8596 0.4711
## Residuals 0.39159 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(Species_counts_cover_rain, aes(x = Site, y = Prop_Immigrant)) +
geom_point()+
labs(x = "", y = "Proportion of immigrant species in seed rain", title = "Seed rain - Aboveground")
# Make a new dataset to log the abundance variables
cover_rain_log <- cover_rain
# Do a log transformation (add the 0.1 to avoid infinite values)
cover_rain_log$Cover <- log(cover_rain_log$Cover+0.1)
cover_rain_log$Number_Seeds <- log(cover_rain_log$Number_Seeds+0.1)
# Get rid of infinite values
cover_rain_log[mapply(is.infinite, cover_rain_log)] <- 0
#cover_rain$Site <- factor(cover_rain$Site, levels = c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))
ggplot(cover_rain_log, aes(x = Cover, y = Number_Seeds))+
geom_point()+
facet_wrap(vars(Site), nrow = 1, strip.position = "top")+
geom_smooth(method = "lm") +
labs(x = "Log(veg. cover)", y = "Log (no. seeds)")+
theme(axis.text = element_text( size = 12),
axis.text.x = element_text( size = 12),
axis.title = element_text( size = 16),
legend.position = "none",
strip.text = element_text(size = 16))
## `geom_smooth()` using formula = 'y ~ x'
ggplot(cover_rain_log, aes(x = Cover, y = Number_Seeds, color = Site))+
geom_point()+
geom_smooth(method = "lm") +
labs(y = "Log (Total Captured Seed Rain )", x = "Log (Total Cover )", title = "Seed rain vs cover")+
theme_classic()
## `geom_smooth()` using formula = 'y ~ x'
# Remove everything that we didn't see both cover and seeds for at a transect
## Removing the "rare" occurrences (cover has to be greater than 1 and number of seeds greater than 10) removes the streaking apparent in the residuals
cover_rain_filtered <- cover_rain %>%
filter(Cover > 0) %>%
filter(Number_Seeds > 0)
Fitting the model
Note I tried other distributional families (gamma, inverse guassian, tweedie, negative binomial) and logging both variables and using a Normal distribution seems to be the best option when looking at the residual vs. fits plot. You do get weird streaks in the residuals though because there were so many instances of us seeing 1 or 2 seeds/cover.
# Fit a linear mixed-model
mod <- lmer(log10(Number_Seeds)~log10(Cover)+(1|Site), data = cover_rain_filtered )
plot(mod)
Histograms of the predictor and response variables after logging show why we get the streaks. They still aren’t particularly Normal after logging.
p <-ggplot(cover_rain_filtered , aes(x=(log10(Number_Seeds)))) +
geom_histogram(color="black", fill="white")+
labs(x = "Log(no. seeds)", y = "Frequency")
p
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
q <-ggplot(cover_rain_filtered , aes(x=log10(Cover))) +
geom_histogram(color="black", fill="white")+
labs(x = "Log(total cover)", y = "Frequency")
q
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Plotting the species and their associated residuals shows which species were over and under represented in the seed rain based on local vegetative cover. It does not show instances where we caught immigrating seeds (no local cover) or species with unsuccessful reproduction (no captured seeds).
resid.cover.rain <- residuals(mod)
SPP6 <- cover_rain_filtered $SPP6
Site <- cover_rain_filtered $Site
resid.species.all.cover.rain <- data.frame(SPP6, resid.cover.rain, Site)
resid.specieslall.cover.rain <- resid.species.all.cover.rain %>%
arrange(desc(resid.cover.rain))
#Of the species captured in both the seed rain and cover, make a plot showing which species were present in the seed rain in greater numbers than expected.
resid.species.all.cover.rain %>%
mutate(SPP6 = fct_reorder(SPP6, resid.cover.rain, .fun='median')) %>%
ggplot(aes(x = SPP6, y = resid.cover.rain, color = Site)) +
labs(x = "Species", y = "Residual", title = "Seed rain vs. Cover")+
geom_point() +
theme_classic() +
theme(axis.text.x = element_text(angle = 70, hjust = 1))+
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1)+
ylim(-5,5)+
theme(axis.text = element_text( size = 12),
axis.text.x = element_text( size = 5),
axis.title = element_text( size = 16),
strip.text = element_text(size = 16)) +
annotate("text", x = 22, y = 3, label = "Overrepresented in seed rain")+annotate("text", x = 75, y = -4, label = "Underrepresented in seed rain")
Here is a plot showing which species were present in the seed rain in greater numbers than expected based on local cover. Bars represent 95% confidence intervals. Absent bars for a species indicates it was only present in the cover and seed rain in one transect.
Remember that this figure only shows species that we were able to captured both in the cover and the seed rain at a transect. Does not show anything we caught only seeds or only cover at.
The figure does line up very well with the results of Maya’s seed predation study. Some of the species favored by rodents were Helianthus mollis, Sporobolus heterolepis, Desmodium canadense, Liatris aspera, and Baptisia bracteata. Based on the figure, these species or congeners also had far fewer seeds that you would expect based on local cover. Other species on this end of the spectrum are also well known to be favored by birds (e.g., American goldfinch and bobwhite). For example, Cirsium altissimum, Silphium integrifolium, Echinacea pallida, Strophostyles leiosperma, and Ambrosia artemisiifolia. Baptisia bracteata and Baptisia alba are also known to suffer heavy pre-dispersal predation from insects at Tucker Prairie (Haddock et al. 1982).
Lysimachia lanceolata is interesting since it has a tight relationship with a rare genera of bees, Macropis. It is unclear how dependent this species is on pollination by Macropis for reproduction.
Size, both plant and seed, are also likely important factors. Based on Moles et al 2004, small plants also tend to produce many small seeds (e.g., Juncus) while large plants (e.g., Baptisia) produce fewer, large seeds with tradeoffs existing between seed size and juvenile survival.
Primary dispersal mode is likely also important here - many species high in cover but low in captured seeds have fleshy fruits. For example, Potentilla simplex, Rhus glabra, Rhus copallinum, Rosa spp., and Fragraria virginica. I believe the only fleshy fruited species we captured was Rubus spp.
Life history strategy is also important, with annuals being grouped on the overrepresented end of the spectrum.
resid.species.all.mean.cover.rain <- resid.species.all.cover.rain %>%
group_by(SPP6) %>%
summarize(mean.resid = mean(resid.cover.rain ),
sd.resid = sd(resid.cover.rain),
n.resid = n()) %>%
mutate(se.resid = sd.resid / sqrt(n.resid),
lower.ci.resid = mean.resid - 1.96*((sd.resid/sqrt(n.resid))),
upper.ci.resid = mean.resid + 1.96*((sd.resid/sqrt(n.resid))))
resid.species.all.mean.cover.rain %>%
mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>%
ggplot(aes(x = SPP6, y = mean.resid)) +
xlab("Species")+
ylab("Residual") +
geom_point() +
geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
theme_classic() +
theme(axis.text.x = element_text(angle = 70, hjust = 1))+
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) +
ylim(-5,5) +
theme(axis.text = element_text( size = 10),
axis.text.x = element_text( size = 5),
axis.title = element_text( size = 16),
strip.text = element_text(size = 16))+
annotate("text", x = 22, y = 1.75, label = "Overrepresented in local seed rain")+
annotate("text", x = 47, y = -1, label = "Underrepresented in local seed rain")
Mass_coverage <- left_join(resid.species.all.mean.cover.rain, traits)
## Joining with `by = join_by(SPP6)`
Mass_coverage <- Mass_coverage %>%
mutate(Seed_Mass = if_else(Mean_1_Seed_Mass_g > 0, 1, 0)) %>%
mutate_all(~replace(., is.na(.), 0))
pane1 <- Mass_coverage %>%
mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>%
ggplot(aes(x = SPP6, y = mean.resid, color = as.factor(Seed_Mass))) +
xlab("Species")+
ylab("Residual") +
geom_point() +
scale_color_manual(values = c("black", "blue"))+
geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
theme_classic() +
theme(axis.text.x = element_text(angle = 70, hjust = 1))+
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) +
ylim(-5,5) +
theme(axis.text = element_text( size = 10),
axis.text.x = element_text( size = 5),
axis.title = element_text( size = 16),
strip.text = element_text(size = 16))+
annotate("text", x = 22, y = 3, label = "Overrepresented in local seed rain")+
annotate("text", x = 75, y = -4, label = "Underrepresented in local seed rain")
pane1
cover_bank_sum <- cover_bank %>%
group_by(SPP6) %>%
summarize(tot.cover = sum(Cover), tot.seeds = sum(Number_Seedlings))
Germinated but didn’t observed in the cover
Almost everything we were able to germinate that was not in the local cover were weedy annual species. Some things like the Cardamine might have been in the cover but have such an early phenology that we missed them. The two tree species might be greenhouse contaminates since the roof was partially open and these trees were right next door.
Observed in the cover but didn’t germinate
Too many to list (112 taxa)
Clearly, only a small fraction of local species survive the transition from cover -> seed rain -> to germinable seed bank.
Few perennial species were present in the germinable seed bank.
## # A tibble: 2,025 × 5
## # Groups: Site, Transect, SPP6 [1,725]
## Site Transect SPP6 Type Abundance
## <chr> <fct> <chr> <chr> <dbl>
## 1 PFCA 1 1 ACAVIR Aboveground 11.5
## 2 PFCA 1 1 ALOCAR Aboveground 1
## 3 PFCA 1 1 AMBART Aboveground 6
## 4 PFCA 1 1 BARVUL Aboveground 3.5
## 5 PFCA 1 1 BROJAP Aboveground 1
## 6 PFCA 1 1 CERSPP Aboveground 1
## 7 PFCA 1 1 CHAFAS Aboveground 388
## 8 PFCA 1 1 CHAFAS Seed_Bank 11
## 9 PFCA 1 1 CONCAN Aboveground 1
## 10 PFCA 1 1 DESILL Aboveground 1
## # … with 2,015 more rows
# Code obtained and modified from Lauren and Anu's seed bank vs. vegetation paper
## Only included Bray-Curtis and Jaccard dissimilarity matrices since the others produced similar results to these two
##### run dissimilarity loops
cover_bank_dis_results <-matrix(nrow=0,ncol=4)
sites <-as.vector(unique(cover_bank_dis$Site))
for(s in 1:length(sites)){
temp <-subset(cover_bank_dis, Site==sites[s])
temp$SPP6 <- factor(as.character(temp$SPP6))
transects <-as.vector(unique(temp$Transect))
#for(b in 1:length(blocks)){
#temp_b<-subset(temp, block==blocks[b])
#plots<-as.vector(unique(temp_b$plot))
for(t in 1:length(transects)){
temp_t <- subset(temp, Transect==transects[t])
#wide form - with a row for above and belowground
plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
#relative abundances
trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"tot"))
names(trans_tot)[1]<-"trt"
trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
names(trans_pa)[1]<-"trt"
distance_bray<-vegdist(trans_tot[,-1], method = "bray")
distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")
new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1])
cover_bank_dis_results <-rbind( cover_bank_dis_results, new.row)
}
print(s/length(sites))
}
## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1
row.names(cover_bank_dis_results )<-NULL
cover_bank_dis_results <- data.frame(cover_bank_dis_results )
names(cover_bank_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard")
##
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = cover_bank_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.220753 -0.058252 0.008957 0.061361 0.165558
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.55757 0.02931 19.023 < 2e-16 ***
## SitePFCA 2 0.26498 0.04145 6.393 2.08e-07 ***
## SitePFCA 3 0.31432 0.04145 7.583 5.73e-09 ***
## SiteTP 0.27515 0.04145 6.638 9.84e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09269 on 36 degrees of freedom
## Multiple R-squared: 0.6679, Adjusted R-squared: 0.6402
## F-statistic: 24.13 on 3 and 36 DF, p-value: 9.766e-09
## Anova Table (Type III tests)
##
## Response: dissimilarity_bray
## Sum Sq Df F value Pr(>F)
## (Intercept) 3.10879 1 361.886 < 2.2e-16 ***
## Site 0.62197 3 24.134 9.766e-09 ***
## Residuals 0.30926 36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call:
## lm(formula = dissimilarity_jaccard ~ Site, data = cover_bank_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.108471 -0.050240 -0.004958 0.054810 0.133643
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.78253 0.02066 37.884 <2e-16 ***
## SitePFCA 2 0.04594 0.02921 1.573 0.1246
## SitePFCA 3 0.06965 0.02921 2.384 0.0225 *
## SiteTP 0.05257 0.02921 1.800 0.0803 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06532 on 36 degrees of freedom
## Multiple R-squared: 0.1474, Adjusted R-squared: 0.07636
## F-statistic: 2.075 on 3 and 36 DF, p-value: 0.1208
## Anova Table (Type III tests)
##
## Response: dissimilarity_jaccard
## Sum Sq Df F value Pr(>F)
## (Intercept) 6.1236 1 1435.2106 <2e-16 ***
## Site 0.0266 3 2.0747 0.1208
## Residuals 0.1536 36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 95% confidence intervals
emmeans::emmeans(Cover_bank.mod.bray, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## PFCA 1 0.558 0.0293 36 0.498 0.617
## PFCA 2 0.823 0.0293 36 0.763 0.882
## PFCA 3 0.872 0.0293 36 0.812 0.931
## TP 0.833 0.0293 36 0.773 0.892
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.558, 0.823, 0.872, 0.833)
lower <- c(0.498, 0.763, 0.812, 0.773)
upper <- c(0.617, 0.882, 0.931, 0.892)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Cover_bank_dissimilarity <- data.frame(site, response, lower, upper)
## 95% confidence intervals
emmeans::emmeans(Cover_bank.mod.jac, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## PFCA 1 0.783 0.0207 36 0.741 0.824
## PFCA 2 0.828 0.0207 36 0.787 0.870
## PFCA 3 0.852 0.0207 36 0.810 0.894
## TP 0.835 0.0207 36 0.793 0.877
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.783, 0.828, 0.852, 0.835)
lower <- c(0.741, .787, 0.810, 0.793)
upper <- c( 0.824, 0.870, 0.894, 0.877)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Cover_bank_jac_dissimilarity <- data.frame(site, response, lower, upper)
# Drop Control for now
cover_bank <- cover_bank[!(cover_bank$Site%in%"CONTROL"),]
# Make a new dataset to log the abundance variables
cover_bank_log <- cover_bank
# Do a log transformation (add the 0.1 to avoid infinite values)
cover_bank_log$Cover <- log(cover_bank_log$Cover+0.1)
cover_bank_log$Number_Seedlings <- log(cover_bank_log$Number_Seedlings+0.1)
# Get rid of infinite values
cover_bank_log[mapply(is.infinite, cover_bank_log)] <- 0
cover_bank_log$Site <- factor(cover_bank_log$Site, levels = c("PFCA 1", "PFCA 2", "PFCA 3", "TP", "CONTROL"), labels = c("Young", "Middle", "Old", "Remnant", "Control"))
ggplot(cover_bank_log, aes(x = Cover, y = Number_Seedlings))+
geom_point()+
facet_wrap(vars(Site), nrow = 1, strip.position = "top")+
geom_smooth(method = "lm") +
labs(x = "Log(veg. cover)", y = "Log (no. seedlings)", title = "Seed bank vs. Cover")+
theme(axis.text = element_text( size = 14),
axis.text.x = element_text( size = 16),
axis.title = element_text( size = 18),
legend.position = "none",
strip.text = element_text(size = 16))
## `geom_smooth()` using formula = 'y ~ x'
ggplot(cover_bank_log, aes(x = Cover, y = Number_Seedlings, color = Site))+
geom_point()+
geom_smooth(method = "lm") +
labs(y = "Log (Total Germinated Seeds )", x = "Log (Total Cover )", title = "Seed bank vs cover")+
theme_classic()
## `geom_smooth()` using formula = 'y ~ x'
Fitting the Model
# Remove everything that we didn't see both cover and seed bank for at a transect
cover_bank_filtered <- cover_bank %>%
filter(Cover > 0) %>%
filter(Number_Seedlings > 0)
#
mod1 <- lmer(log(Number_Seedlings)~log(Cover)+(1|Site), data = cover_bank_filtered )
plot(mod1)
Again we get the streaks in the residual plots.
p1<-ggplot(cover_bank_filtered , aes(x=log(Number_Seedlings))) +
geom_histogram(color="black", fill="white")
p1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
q1 <-ggplot(cover_bank_filtered , aes(x=log(Cover))) +
geom_histogram(color="black", fill="white")
q1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
resid.cover.bank <- residuals(mod1)
#resid <- filtered.mod$residuals
SPP6 <- cover_bank_filtered $SPP6
Site <- cover_bank_filtered $Site
resid.species.all.cover.bank <- data.frame(SPP6, resid.cover.bank, Site)
resid.specieslall.cover.bank <- resid.species.all.cover.bank %>%
arrange(desc(resid.cover.bank))
#Of the species captured in both the seed bank and cover, make a plot showing which species were present in the seed bank in greater numbers than expected.
resid.species.all.cover.bank %>%
mutate(SPP6 = fct_reorder(SPP6, resid.cover.bank, .fun='median')) %>%
ggplot(aes(x = SPP6, y = resid.cover.bank, color = Site)) +
labs(x = "Species", y = "Residual", title = "Seed Bank vs. Cover")+
geom_point() +
theme_classic() +
theme(axis.text.x = element_text(angle = 70, hjust = 1))+
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) +
theme(axis.text = element_text( size = 10),
axis.text.x = element_text( size = 8),
axis.title = element_text( size = 16),
strip.text = element_text(size = 16)) +
ylim(-4,4)+
annotate("text", x = 12, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 34, y = -4, label = "Underrepresented in seed bank")
Looking at this residual plot, disturbance tolerant annual species tend to be over represented in the seed bank compared to the local cover. These species are most likely enriched in seed banks overtime and are persistent. For example, Digitaria ischaemum, Setaria faberi, Plantago virginica, and Juncus spp. Some over representation might be because of superb dispersal ability (Eupatorium serotinum and Agrostis hyemalis).
Life history is clearly important, with annuals being more present than perennials.
The seed bank and cover had the least amount of overlap with shared species.
resid.species.all.mean.cover.bank <- resid.species.all.cover.bank %>%
group_by(SPP6) %>%
summarize(mean.resid = mean(resid.cover.bank),
sd.resid = sd(resid.cover.bank),
n.resid = n()) %>%
mutate(se.resid = sd.resid / sqrt(n.resid),
lower.ci.resid = mean.resid - 1.96*((sd.resid/sqrt(n.resid))),
upper.ci.resid = mean.resid + 1.96*((sd.resid/sqrt(n.resid))))
resid.species.all.mean.cover.bank %>%
mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>%
ggplot(aes(x = SPP6, y = mean.resid)) +
labs(x = "Species", y = "Residual", title = "Seed Bank vs. Cover")+
geom_point() +
geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
theme_classic() +
theme(axis.text.x = element_text(angle = 70, hjust = 1))+
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) + theme(axis.text = element_text( size = 12),
axis.text.x = element_text( size = 8),
axis.title = element_text( size = 16),
strip.text = element_text(size = 16))+
ylim(-4,4)+
annotate("text", x = 10, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 34, y = -4, label = "Underrepresented in seed bank")
Mass_coverage_cover.bank <- left_join(resid.species.all.mean.cover.bank, traits)
## Joining with `by = join_by(SPP6)`
Mass_coverage_cover.bank <- Mass_coverage_cover.bank %>%
mutate(Seed_Mass = if_else(Mean_1_Seed_Mass_g > 0, 1, 0)) %>%
mutate_all(~replace(., is.na(.), 0))
pane2 <- Mass_coverage_cover.bank %>%
mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>%
ggplot(aes(x = SPP6, y = mean.resid, color = as.factor(Seed_Mass))) +
xlab("Species")+
ylab("Residual") +
geom_point() +
scale_color_manual(values = c("black", "blue"))+
geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
theme_classic() +
theme(axis.text.x = element_text(angle = 70, hjust = 1))+
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) +
ylim(-5,5) +
theme(axis.text = element_text( size = 10),
axis.text.x = element_text( size = 5),
axis.title = element_text( size = 16),
strip.text = element_text(size = 16))+
annotate("text", x = 16, y = 3, label = "Overrepresented in local seed bank")+
annotate("text", x = 32, y = -4, label = "Underrepresented in local seed bank")
pane2
rain_bank_sum <- rain_bank %>%
group_by(SPP6) %>%
summarize(tot.seeds = sum(Number_Seeds), tot.seedlings = sum(Number_Seedlings))
Germinated but didn’t catch in the seed rain
Mostly ruderal annual species.
Caught in the seed rain but didn’t germinate
Too many to list (61 taxa)
## # A tibble: 1,859 × 5
## # Groups: Site, Transect, SPP6 [1,500]
## Site Transect SPP6 Type Abundance
## <chr> <fct> <chr> <chr> <dbl>
## 1 PFCA 1 1 ACAVIR Seed_Rain 22
## 2 PFCA 1 1 AMBART Seed_Rain 6
## 3 PFCA 1 1 ANAMIN Seed_Rain 11
## 4 PFCA 1 1 BOLAST Seed_Rain 1
## 5 PFCA 1 1 CERSPP Seed_Rain 14
## 6 PFCA 1 1 CHAFAS Seed_Bank 11
## 7 PFCA 1 1 CHAFAS Seed_Rain 276
## 8 PFCA 1 1 CHEALB Seed_Rain 1
## 9 PFCA 1 1 CONCAN Seed_Rain 5
## 10 PFCA 1 1 DIGISC Seed_Bank 1
## # … with 1,849 more rows
# Code obtained and modified from Lauren and Anu's seed bank vs. vegetation paper
## Only included Bray-Curtis and Jaccard dissimilarity matrices since the others produced similar results to these two
##### run dissimilarity loops
rain_bank_dis_results <-matrix(nrow=0,ncol=4)
sites <-as.vector(unique(rain_bank_dis$Site))
for(s in 1:length(sites)){
temp <-subset(rain_bank_dis, Site==sites[s])
temp$SPP6 <- factor(as.character(temp$SPP6))
transects <-as.vector(unique(temp$Transect))
for(t in 1:length(transects)){
temp_t <- subset(temp, Transect==transects[t])
#wide form - with a row for above and belowground
plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
#relative abundances
trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"tot"))
names(trans_tot)[1]<-"trt"
trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
names(trans_pa)[1]<-"trt"
distance_bray<-vegdist(trans_tot[,-1], method = "bray")
distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")
new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1])
rain_bank_dis_results <-rbind(rain_bank_dis_results, new.row)
}
print(s/length(sites))
}
## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1
row.names(rain_bank_dis_results)<-NULL
rain_bank_dis_results <- data.frame(rain_bank_dis_results)
names(rain_bank_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard")
##
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = rain_bank_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30162 -0.06832 -0.00551 0.07355 0.26438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.55122 0.04061 13.575 9.94e-16 ***
## SitePFCA 2 0.24085 0.05743 4.194 0.00017 ***
## SitePFCA 3 0.18002 0.05743 3.135 0.00342 **
## SiteTP 0.13699 0.05743 2.386 0.02244 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1284 on 36 degrees of freedom
## Multiple R-squared: 0.3458, Adjusted R-squared: 0.2913
## F-statistic: 6.344 on 3 and 36 DF, p-value: 0.001449
## Anova Table (Type III tests)
##
## Response: dissimilarity_bray
## Sum Sq Df F value Pr(>F)
## (Intercept) 3.03840 1 184.2741 9.941e-16 ***
## Site 0.31379 3 6.3436 0.001449 **
## Residuals 0.59359 36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call:
## lm(formula = dissimilarity_jaccard ~ Site, data = rain_bank_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.135047 -0.042571 0.002154 0.030196 0.128868
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.72479 0.02097 34.558 < 2e-16 ***
## SitePFCA 2 0.02556 0.02966 0.862 0.39458
## SitePFCA 3 0.08522 0.02966 2.873 0.00677 **
## SiteTP 0.04026 0.02966 1.357 0.18317
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06632 on 36 degrees of freedom
## Multiple R-squared: 0.1949, Adjusted R-squared: 0.1278
## F-statistic: 2.905 on 3 and 36 DF, p-value: 0.04791
## Anova Table (Type III tests)
##
## Response: dissimilarity_jaccard
## Sum Sq Df F value Pr(>F)
## (Intercept) 5.2532 1 1194.2473 < 2e-16 ***
## Site 0.0383 3 2.9051 0.04791 *
## Residuals 0.1584 36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 95% confidence intervals
emmeans::emmeans(Rain_bank.mod.bray, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## PFCA 1 0.551 0.0406 36 0.469 0.634
## PFCA 2 0.792 0.0406 36 0.710 0.874
## PFCA 3 0.731 0.0406 36 0.649 0.814
## TP 0.688 0.0406 36 0.606 0.771
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.551, 0.792, 0.731, 0.688)
lower <- c(0.469, 0.710, 0.649, 0.606)
upper <- c(0.634, 0.874, 0.814, 0.771)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Rain_bank_dissimilarity <- data.frame(site, response, lower, upper)
## 95% confidence intervals
emmeans::emmeans(Rain_bank.mod.jac, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## PFCA 1 0.725 0.021 36 0.682 0.767
## PFCA 2 0.750 0.021 36 0.708 0.793
## PFCA 3 0.810 0.021 36 0.767 0.853
## TP 0.765 0.021 36 0.723 0.808
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.725, 0.750, 0.810, 0.765)
lower <- c(0.682, 0.708, 0.767, 0.723)
upper <- c(0.767, 0.793, 0.853, 0.808)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Rain_bank_jac_dissimilarity <- data.frame(site, response, lower, upper)
rain_bank_dis_wide <- rain_bank_dis %>%
spread(SPP6, Abundance) %>%
mutate_all(~replace(., is.na(.), 0))
## `mutate_all()` ignored the following grouping variables:
## • Columns `Site`, `Transect`
## ℹ Use `mutate_at(df, vars(-group_cols()), myoperation)` to silence the message.
rain_bank_dis_long <- rain_bank_dis_wide %>%
gather(key = "SPP6", value = "Abundance", 4:144) %>%
spread(key = Type, Abundance)
# Get rid of columns that both have zeros
rain_bank_dis_long_unique <- rain_bank_dis_long[!(rain_bank_dis_long$Seed_Bank == 0 & rain_bank_dis_long$Seed_Rain == 0),]
# 0 = no seeds/seedlings and 1 = seeds/seedlings
rain_bank_dis_long_unique <- rain_bank_dis_long_unique %>%
mutate(Seed_Bank = if_else(Seed_Bank > 0, 1, 0)) %>%
mutate(Seed_Rain = if_else(Seed_Rain > 0, 1, 0))
# Make a new column that adds both the seed rain and seed bank presence/absence columns. Then remove rows that = 2 (shares both species)
rain_bank_dis_long_unique$shared <- rain_bank_dis_long_unique$Seed_Bank + rain_bank_dis_long_unique$Seed_Rain
rain_bank_dis_long_unique2 <- rain_bank_dis_long_unique[!(rain_bank_dis_long_unique$shared == 2),]
rain_bank_dis_long_unique2 <- rain_bank_dis_long_unique2[, -6]
# Count unique taxa per site and transect for the seed bank compared to the seed rain
Unique_species_bank_com_rain <- rain_bank_dis_long_unique2[, -5]
Unique_species_bank_com_rain <- Unique_species_bank_com_rain %>%
filter(Seed_Bank != 0 ) %>% group_by(Site, Transect) %>%
summarize(New_Bank = length(unique(SPP6)))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
# Count unique taxa per site and transect for the seed rain compared to the seed bank
Unique_species_rain_com_bank <- rain_bank_dis_long_unique2[, -4]
Unique_species_rain_com_bank <- Unique_species_rain_com_bank %>%
filter(Seed_Rain != 0 ) %>%
group_by(Site, Transect) %>%
summarize(New_Rain = length(unique(SPP6)))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
# Count the number of shared taxa between communities
rain_bank_dis_long_shared <- rain_bank_dis_long_unique %>%
filter(shared == 2) %>%
group_by(Site, Transect) %>%
summarize(Shared = length(unique(SPP6)))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
## Merge them together
Species_counts_rain_bank <- full_join(Unique_species_bank_com_rain, Unique_species_rain_com_bank) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:
## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 3 remaining warnings.
Species_counts_rain_bank <- full_join(Species_counts_rain_bank, rain_bank_dis_long_shared) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:
## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 3 remaining warnings.
Species_counts_rain_bank$Tot_Richness <- Species_counts_rain_bank$New_Bank + Species_counts_rain_bank$New_Rain + Species_counts_rain_bank$Shared
# Count the number of species in the seed rain for each transect
Species_counts_rain_bank$Tot_Rain <- Species_counts_rain_bank$New_Rain + Species_counts_rain_bank$Shared
# Count the number of species in the seed bank for each transect
Species_counts_rain_bank$Tot_Bank <- Species_counts_rain_bank$New_Bank + Species_counts_rain_bank$Shared
## Immigrants
Species_counts_rain_bank$Prop_Immigrant <- Species_counts_rain_bank$New_Bank / Species_counts_rain_bank$Tot_Rain
ggplot(Species_counts_rain_bank, aes(x = Site, y = Prop_Immigrant)) +
geom_point()
# Drop the control for now
rain_bank <- rain_bank[!(rain_bank$Site%in%"CONTROL"),]
# Make a new dataset to log the abundance variables
rain_bank_log <- rain_bank
# Do a log transformation (add the 0.1 to avoid infinite values)
rain_bank_log$Number_Seedlings <- log(rain_bank_log$Number_Seedlings+0.1)
rain_bank_log$Number_Seeds <- log(rain_bank_log$Number_Seeds+0.1)
# Get rid of infinite values
rain_bank_log[mapply(is.infinite, rain_bank_log)] <- 0
rain_bank_log$Site <- factor(rain_bank_log$Site, levels = c("PFCA 1", "PFCA 2", "PFCA 3", "TP", "CONTROL"), labels = c("Young", "Middle", "Old", "Remnant", "Control"))
ggplot(rain_bank_log, aes(x = Number_Seeds, y = Number_Seedlings))+
geom_point()+
facet_wrap(vars(Site), nrow = 1, strip.position = "top")+
geom_smooth(method = "lm") +
labs(x = "Log(no. seeds)", y = "Log (no. seedlings)")+
theme(axis.text = element_text( size = 14),
axis.text.x = element_text( size = 20),
axis.title = element_text( size = 18),
legend.position = "none",
strip.text = element_text(size = 16))
## `geom_smooth()` using formula = 'y ~ x'
ggplot(rain_bank_log, aes(x = Number_Seeds, y = Number_Seedlings, color = Site))+
geom_point()+
geom_smooth(method = "lm") +
labs(y = "Log (Total Captured Seed Rain )", x = "Log (Total Cover )", title = "Seed rain vs Seed Bank")+
theme_classic()
## `geom_smooth()` using formula = 'y ~ x'
# Remove everything that we didn't see both the seed rain and bank at a transect
rain_bank_filtered <- rain_bank %>%
filter(Number_Seedlings> 0) %>%
filter(Number_Seeds > 0)
#
mod3 <- lmer(log(Number_Seedlings)~log(Number_Seeds)+(1|Site), data = rain_bank_filtered)
plot(mod3)
Again we get the streaks in the residual vs. fitted values plot.
p2<-ggplot(rain_bank_filtered , aes(x=log(Number_Seedlings))) +
geom_histogram(color="black", fill="white")
p2
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
q2 <-ggplot(rain_bank_filtered , aes(x=log(Number_Seeds))) +
geom_histogram(color="black", fill="white")
q2
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
resid.rain.bank<- residuals(mod3)
SPP6 <- rain_bank_filtered$SPP6
Site <- rain_bank_filtered$Site
resid.species.all.rain.bank <- data.frame(SPP6, resid.rain.bank, Site)
resid.specieslall.rain.bank <- resid.species.all.rain.bank %>%
arrange(desc(resid.rain.bank))
#Of the species captured in both the seed rain and bank, make a plot showing which species were present in the seed bank in greater numbers than expected.
resid.species.all.rain.bank %>%
mutate(SPP6 = fct_reorder(SPP6, resid.rain.bank, .fun='median')) %>%
ggplot(aes(x = SPP6, y = resid.rain.bank, color = Site)) +
labs(x = "Species", y = "Residual", title = "Seed Rain vs. Seed Bank") +
geom_point() +
theme_classic() +
theme(axis.text.x = element_text(angle = 70, hjust = 1))+
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1)+
theme(axis.text = element_text( size = 12),
axis.text.x = element_text( size = 8),
axis.title = element_text( size = 16),
strip.text = element_text(size = 16))+
ylim(-4,4)+
annotate("text", x = 12, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 40, y = -4, label = "Underrepresented in seed bank")
resid.species.all.mean.rain.bank <- resid.species.all.rain.bank %>%
group_by(SPP6) %>%
summarize(mean.resid = mean(resid.rain.bank),
sd.resid = sd(resid.rain.bank),
n.resid = n()) %>%
mutate(se.resid = sd.resid / sqrt(n.resid),
lower.ci.resid = mean.resid - 1.96*((sd.resid/sqrt(n.resid))),
upper.ci.resid = mean.resid + 1.96*((sd.resid/sqrt(n.resid))))
Interesting that some species underrepresented in the seed rain based on local cover are present in germinable seed bank (LYSLAN, CIRALT, STRLEI, EUTGYM, RATPIN). Does this suggest high mortality when transitionig into the seed bank and low mortality once in the seed bank?
Again, most species overrepresented in the seed bank based on local seed input are weedy and annual species. Most perennials fall on the underrepresented side of the spectrum.
Many of the native ones also seem to be species that do well in restorations (Lespedeza capitata, Ratibida pinnata, Andropogon gerardii, Sorghastrum nutans, Coreopsis lanceolata, Tradescantia ohiensis, Penstemon digitalis, Schizachyrium scoparium, Solidago altissima, Symphytrichum spp., etc.). It’s probably no coincidence that these species that recruit from seed rain also do well in restorations that are broadcast seeded. I wonder if restoration practices are accidentally selecting for species that recruit well from seed (limited dispersal opportunities and reinforcement from what initially recruits).
I think the most surprising species that germinated is Lysimachia lanceolata, which is not present in the PFCA restorations.
resid.species.all.mean.rain.bank %>%
mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>%
ggplot(aes(x = SPP6, y = mean.resid)) +
labs(x = "Species", y = "Residual", title = "Seed Rain vs. Seed Bank") +
geom_point() +
geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
theme_classic() +
theme(axis.text.x = element_text(angle = 70, hjust = 1))+
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) + theme(axis.text = element_text( size = 12),
axis.text.x = element_text( size = 8),
axis.title = element_text( size = 16),
strip.text = element_text(size = 16))+
ylim(-4,4)+
annotate("text", x = 12, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 40, y = -4, label = "Underrepresented in seed bank")
resid.species.all.mean.rain.bank %>%
mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>%
ggplot(aes(x = SPP6, y = mean.resid)) +
labs(x = "Species", y = "Residual", title = "Seed Bank vs. Cover")+
geom_point() +
geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
theme_classic() +
theme(axis.text.x = element_text(angle = 70, hjust = 1))+
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) + theme(axis.text = element_text( size = 12),
axis.text.x = element_text( size = 8),
axis.title = element_text( size = 16),
strip.text = element_text(size = 16))+
ylim(-4,4)+
annotate("text", x = 10, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 34, y = -4, label = "Underrepresented in seed bank")
Mass_coverage_rain.bank <- left_join(resid.species.all.mean.rain.bank, traits)
## Joining with `by = join_by(SPP6)`
Mass_coverage_rain.bank <- Mass_coverage_cover.bank %>%
mutate(Seed_Mass = if_else(Mean_1_Seed_Mass_g > 0, 1, 0)) %>%
mutate_all(~replace(., is.na(.), 0))
pane3 <- Mass_coverage_rain.bank %>%
mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>%
ggplot(aes(x = SPP6, y = mean.resid, color = as.factor(Seed_Mass))) +
xlab("Species")+
ylab("Residual") +
geom_point() +
scale_color_manual(values = c("black", "blue"))+
geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
theme_classic() +
theme(axis.text.x = element_text(angle = 70, hjust = 1))+
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) +
ylim(-5,5) +
theme(axis.text = element_text( size = 10),
axis.text.x = element_text( size = 5),
axis.title = element_text( size = 16),
strip.text = element_text(size = 16))+
annotate("text", x = 16, y = 3, label = "Overrepresented in local seed rain")+
annotate("text", x = 32, y = -4, label = "Underrepresented in local seed rain")
pane3